module polyroot
    implicit none
    integer,parameter :: N=3
    real,parameter :: poly(N+1)=(/-3.0,11.0,-6.0,1.0/)
contains
    function f(x)
        implicit none
        real :: f,x,t
        integer :: i
        t=0.0
        do i=1,N+1
            t=t*x+poly(N-i+2)
        end do
        f=t
    end function
    subroutine root(r)
        implicit none
        real,allocatable :: r(:)
        real :: p(N+1),limit,M
        allocate(r(N))
        p=poly/poly(N+1)
        limit=sum(abs(p(:N)))
        if(limit < 1.0)then
            limit=1.0
        end if
        call ddhrt(-limit,limit,limit/50.0,1d-5,r,N,M,f)
    end subroutine
end module
program main
    use polyroot
    implicit none
    real,allocatable :: r(:)
    call root(r)
    write(*,*)r
    deallocate(r)
end program

SUBROUTINE DDHRT( A , B, H , EPS, X, N , M , F )
    DIMENSION X( N )
    M= 0
    Z= A
    Y = F(Z)
10  IF ( ( Z .GT. B+ H / 2.0) .OR. ( M .EQ. N ) )RETURN
    IF( ABS( Y ) .LT. EPS) THEN
        M= M+ 1
        X( M ) = Z
        Z= Z+ H / 2.0
        Y = F ( Z)
        GOTO 10
    END IF
    Z1= Z + H
    Y1= F (Z1)
    IF ( ABS( Y1) .LT. EPS ) THEN
        M= M+ 1
        X( M ) = Z1
        Z= Z1+ H / 2.0
        Y = F ( Z)
        GOTO 10
    END IF
    IF ( Y * Y1.GT. 0.0) THEN
        Y= Y1
        Z= Z1
        GOTO 10
    END IF
20  IF ( ABS( Z1- Z ) .LT. EPS) THEN
        M= M+ 1
        X( M ) = ( Z1+ Z) / 2.0
        Z = Z1+ H / 2.0
        Y = F ( Z)
        GOTO 10
    END IF
    Z0= ( Z1+ Z ) / 2.0
    Y0= F ( Z0)
    IF ( ABS( Y0) .LT. EPS ) THEN
        M= M+ 1
        X( M ) = Z0
        Z= Z0+ H / 2.0
        Y = F ( Z)
        GOTO 10
    END IF
    IF ( Y * Y0.LT. 0.0) THEN
        Z1= Z0
        Y1= Y0
    ELSE
        Z= Z0
        Y= Y0
    END IF
    GOTO 20
END
